import numpy as np
from scipy.constants import codata
from random import uniform
class OneD_Ising:
	def __init__(self, n = 0, J = 1., T = 0., B = 0.):
		self.up = 1.
		self.down = -1.
		self.J = J
		self.n = n
		self.B = B
		self.chain = np.zeros(self.n)
		self.currH = 0.
		self.sumH = 0.
		self.beta = 1/float(T)

	def Thermalize(self):
		for i in range(self.n):
			self.chain[i] = self.up if int(round(uniform(0, 1))) else self.down
		self.currH = self.GetEnergy()

	def GetEnergy(self):
		H = 0.
		for i in range(self.n):
			# periodic boundary conditions
			H += -self.J*self.chain[i]*self.chain[i-1] - self.B*self.chain[i]
		return H

	def FlipSpin(self, n):
		self.chain[n] = self.up if self.chain[n] == self.down else self.down

	def LocalUpdate(self):
		tmpH = 0.
		n = int(round(uniform(0, 1)*(self.n-1)))
		self.FlipSpin(n)
		tmpH = self.GetEnergy()
		if tmpH > self.currH:
			if uniform(0, 1) > np.e**((self.currH-tmpH)*self.beta):
				self.FlipSpin(n) #undo
			else:
				self.currH = tmpH
		else:
			self.currH = tmpH
		self.sumH += self.currH

	def Simulate(self, N):
		self.currH = 0.
		self.sumH = 0.
		self.Thermalize()
		for i in range(N):
			self.LocalUpdate()
		print self.chain
		print self.currH
		print self.sumH/float(N)


T = [float(i+1)*10. for i in range(10)]
N = 10000
Ts = []
Hs = []
for i in T:
	t = OneD_Ising(20, 1., i, 1.)
	t.Simulate(N)
	Ts.append(i)
	Hs.append(t.sumH/float(N))

import Gnuplot
g = Gnuplot.Gnuplot()
g.plot(Hs)
raw_input()
